DOE/ER/41132-100-INT00 
UCRL-JC-142857 



Observing Non-Gaussian Sources in Heavy-Ion Reactions 

D.A. Brown 1 . P. Danielewicz 2 
1 University of Washington, Seattle, Washington 98195 and 
Lawrence Livermore National Laboratory, Livermore California 94551 
2 Michigan State University, East Lansing, Michigan 48824 
(February 8, 2008) 



o 
o 



(N 
> 
oo 
o 



o 

O 



o 

=5 



X 



We examine the possibility of extracting non-Gaussian 
sources from two-particle correlations in heavy-ion reactions. 
Non-Gaussian sources have been predicted in a variety of 
model calculations and may have been seen in various like- 
meson pair correlations. As a tool for this investigation, we 
have developed an improved imaging method that relies on a 
Basis spline expansion of the source functions with an im- 
proved implementation of constraints. We examine under 
what conditions this improved method can distinguish be- 
tween Gaussian and non-Gaussian sources. Finally, we in- 
vestigate pion, kaon, and proton sources from the p-Pb re- 
action at 450 GeV/nucleon and from the S-Pb reaction at 
200 GeV/nucleon studied by the NA44 experiment. Both the 
pion and kaon sources from the S-Pb correlations seem to ex- 
hibit a Gaussian core with an extended, non-Gaussian halo. 
We also find evidence for a scaling of the source widths with 
particle mass in the sources from the p-Pb reaction. 



I. INTRODUCTION 

Two-particle correlations have proven to be an impor- 
tant tool for experimentally accessing the space-time ex- 
tent of heavy-ion collisions. For like-meson pair corre- 
lations (e.g. 7r's and K's), the correlation is dominated 
by the so-called Hanbury-Brown/Twiss (HBT) effect (in 
other words, Bosc-Einstcin symmctrization of the meson- 
pair wavefunction) and the Coulomb corrected corre- 
lations are usually adequately parameterized by Gaus- 
sians Since meson final state interactions (FSI) 

can usually be neglected, the Coulomb corrected correla- 
tion function becomes very nearly the Fourier transform 
of a source function. Thus, a Gaussian correlation corre- 
sponds to a Gaussian source function. In general, there 
is no reason to expect the source to be Gaussian. In fact, 
non-Gaussian sources may already have been observed in 
data @4T|. 

There are several reasons to expect non-Gaussian 
sources: contributions from resonance decays should 
lead to an exponential halo Jl3|-p^|, effects of space- 



momentum correlations (caused by either flow [Q,^6[ or 
string fragmentation jl?]]) should lead to a focusing of 
the source jl6j , and even simple geometry should lead to 
non-Gaussian sources. Experimentally distinguishing be- 
tween Gaussian and non-Gaussian sources is difficult and 
may be complicated by FSI within the pair. Recently it 
was realized that, by applying imaging techniques to the 
correlation data, we may extract the two-particle source 
function directly fUJjJ. The imaging has two main ad- 
vantages over the traditional HBT approach: it is model 
independent, meaning that it may reveal non-Gaussian 
features in the source, and it can clearly separate the 
effects of the FSI and symmctrization from effects due 
to the source itself. This last point requires elaboration: 
imaging extracts source functions that may be directly 
compared using correlations that cannot easily be com- 
pared when they arise from completely different parti- 
cles. It is this feature that allows us to compare proton, 
kaon, and pion sources from p-Pb and S-Pb reactions 
from NA44. Indeed, a direct comparison of the proton 
and kaon sources from the p-Pb reaction suggests a sim- 
ple scaling of the source widths that one should expect 
based on Lund-type string phenomenology in a fragment- 
ing string. 

Extracting the source function, Sp(r'), begins by not- 
ing that Sp(r') is related to the experimentally measured 
two-particle correlation, Cp(q'), through a simple linear 
integral equation 0,|l8| : 



Cp(q')-1= / dr'K{d',v')S P (T'). 



(1) 



Thus, "imaging the source" means somehow inverting 
this equation. Here primes denote quantities in the pair 
center of mass (CM) frame. Although for imaging pur- 
poses it is simplest to write Eq. ([!]) in the pair CM, 
Eq. (|l]) may be written in any frame as lZp(q') is a 
Lorentz invariant observable. In (|I|), P = pi + P2 is 
the total momentum of the pair in the lab frame. The P 
subscript indicates the boost from the lab to the pair CM 
frame (P/Pq is the boost velocity between the frames). 
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The kernel of Eq. (m) is 



if(q',r') = |<7 ) (r')| 2 -l. 



(2) 



The wavefunction, describes the propagation of the 

pair from a relative separation of r' in the pair CM to the 
detector with relative momentum q' = |(pi — p' 2 ). The 
source function itself is the quasi-probability of emitting 
the pair a distance of r' apart, in the CM frame. We 
write the source as a convolution of Wigner functions, 
D(r,t,P/2): 



S P (r') = 

J dt'J d 3 RdT D(R + r/2,T + t/2,P/2) 
x D(R-r/2,T-i/2,P/2), 



(3) 



where the variables in the lab frame are understood as 
functions of the variables in the pair CM frame. Here the 
Wigner functions are normalized particle emission rates 



D(r,t,p) 



Ed 7 N 



Ed 3 N 



d 3 r dt d 3 p I d 3 p 



(4) 



and may be computed directly from a transport model as 
discussed in [pTln2 16 1. Due to the time integral in (0), 



we cannot distinguish whether a given r' is associated 
with a time separation or a spatial separation. 

Inverting Eq. ([!]) is generally an ill-posed problem. 
This means that small fluctuations in the data, even 
if well within statistical or systematic errors, can lead 
to large changes in the imaged source function. 111- 
posedness stems from experimental factors (e.g. limited 
statistics, finite sized momentum bins, etc.) and the in- 
trinsic resolution of the kernel in Eq. (0). In other fields, 
this stability problem is attacked using a variety of tac- 
tics including forcing the source function to obey known 
constraints or choosing a representation of the problem 
in which the kernel's resolution may be optimized. Both 
of these techniques were exploited in Ref. j^] . While the 
imaging in Ref. jl2| was successful, the restored sources 
were represented in a basis that does not exhibit the con- 
tinuity that we expect to see in the source. In this paper, 
we report a dramatic improvement of the imaging by us- 
ing a representation of the source in which wc have direct 
control over the continuity of the source. Our choice of 
representation still allows us to utilize constraints and to 
optimize the resolution of the kernel. 

This paper is organized as follows. First, we will set up 
the problem of inverting angle-averaged correlations (i.e. 
expressed in terms of qt nv = ^/q 2 — 9o ) ana - outline the 
improved imaging method. The details of the imaging 
method and our representation of the source are con- 
tained in the appendices. Next, we apply the imaging 
method to correlations corresponding to Gaussian and 
non-Gaussian sources. This will orient us to some of the 
issues we will face when examining real data. Finally, we 



will confront like-pion, likc-kaon, and two-proton correla- 
tion data from S-Pb collisions at 200 GeV/nucleon from 
NA44 and p+Pb collisions at 450 GeV/nucleon. 



II. STATEMENT OF THE PROBLEM 

In this section, we will set up the imaging problem. 
To simplify our discussion, we will consider only angle- 
averaged correlations and sources. First, we will outline 
the one-dimensional imaging problem and mention some 
of our expectations based on experience with Fourier 
transforms. Second, wc will outline how wc utilize the 
Basis spline representation. Finally, we will describe our 
solution using a Bayesian approach to imaging. The de- 
tails of the Basis spline basis and the imaging itself are 
included in the appendices. 



A. The one-dimensional imaging problem 

The source function and correlation in Eq. ([!]) both 
may be expanded in spherical harmonics and the re- 
lations between the angular coefficients are listed in 
Ref. Q. With this expansion, we may image the individ- 
ual components of the correlation function and compile a 
full three-dimensional imaged source. When doing this, 
one must take care in interpreting the results: imaging 
is formulated in the pair CM frame as opposed to frame 
in which we have more intuition, e.g. the lab frame. In 
this letter, we work with only the first term in the spher- 
ical expansions, i.e. the angle-averaged source and cor- 
relation. The drawback of performing a one-dimensional 
analysis is that the angular information is lost and the 
resulting source function is even more difficult to inter- 
pret. 

The angle averaged version of Eq. (fy) is 



K(q) = C{q) - 1 = 4tt J dr r 2 K{q,r)S{r) 



(5) 



Here q = |q'|. For like pairs in the pair CM frame, q' Q = 0. 
This implies that |q'| = q inv = \/q 2 — 9q ■ 

In Eq. (H), the kernel is simply the kernel in Eq. (||), 
but averaged over the angle between q and r: 



1 f 1 



(6) 



For identical spin-zero bosons with no FSI, this kernel is 
K(q,r) = sin (2qr)/2qr, (7) 
while with FSI it is, 



K(q,r)= £ 



\9i(r)\ 2 
{2£+l) 



(8) 
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Here I is the orbital angular momentum quantum num- 
ber. Finally, for protons, the spin-averaged kernel is 



K( q ,r) = l^2(2j + l)(g%(r)Y-l 



jsW 



(9) 



Here t and £' are both orbital angular momentum quan- 
tum numbers and j and s are the total angular momen- 
tum and spin quantum numbers. In the last two cases, g 
is the relative final-state radial wavefunction. For uncor- 
rected meson data, g is the solution of the Klein-Gordon 
equation including the Coulomb potential. For protons, 
g is the solution of the Schrodingcr equation using the 
Coulomb potential and REID93 |nj nucleon-nucleon po- 
tential. 

Given that the identical particle kernels in Eq. (|l|) 
or (||) are Fourier transform kernels at large distances, 
we expect our transforms to behave like Fourier trans- 
forms. If Eq. (5) were a Fourier transform, then by dis- 
cretizing Eq. (5), we would be converting the imaging 
problem into a finite Fourier transform. In this case, the 
Sampling Theorem tells us how the sizes of the bins and 
the numbers of bins in the appropriate Fourier spaces are 
related: 

heir , . heir 

Aq = (10) 



and 



Qmax 1 max 

where q ma x = NAq, r max = NAr and N is the number 
of bins in both r and q space. 

Using these relations, we may get a feeling for how 
structure in the data effect the imaged source. For exam- 
ple, the low-g structure in the data sets the large length 
scale behavior of the source. Conversely, the high-g por- 
tion of the data sets the short length scale behavior of the 
source and therefore sets the size of the smallest features 
features we could hope to resolve in the source. For ex- 
ample, if the correlation dies off around a q f» 80 MeV/c, 
then we should not expect to resolve structure smaller 
than Ar ss 8 fm. Owing to the fact that our kernel is not 
a trigonometric function in general, these estimates are 
qualitative at best. Nevertheless, we will often appeal to 
Fourier theory to explanations of some of the effects that 
we see while imaging. 



B. The representation of the problem 

In our calculations, we expand the imaged source in a 
function basis 



(11) 



and, in this basis, the error on the source is 
AS(r) 



N M 



]T A^SijB^Bjir). (12) 



Here A 2 S is the covariance matrix of the source coeffi- 
cients. Once we average the kernel over momentum bins 
to account for the experimental binning, our inversion 
problem reduces to the following matrix equation: 



N M 
3 = 1 

where the kernel matrix is 

^ n pqi+Aq/2 roc 

A" / d 1 / drr^Kiq^B^r), 

^1 Jq t -Aq/2 JO 



(13) 



(14) 



Here Aq is the momentum bin size. Our source vector 
is made of the coefficients Sj of the basis function rep- 
resentation of the source and our data vector is made of 
the correlation values IZi. 

The function basis which we use to represent our source 
function must have several properties: 1) it must be effi- 
cient, i.e. requiring few coefficients to represent a realistic 
source, 2) it should be continuous, or at least have con- 
tinuity as an option, and 3) it should have an adjustable 
parameter that we might use to optimize the resolution 
in a manner analogous to Ref. |12| . One obvious possi- 
bility is to either use a Laguerre expansion (2^] (so that 
the first term is an exponential fitted to the source) or an 
Edgeworth expansion [^o|,0 (so that the leading term is 
a Gaussian fitted to the source). The downside of either 
of these choices is that it difficult to adjust the terms in 
one's expansion to maximize the resolution of the inver- 
sion. Furthermore, one could argue that if one picks one 
of these bases and keeps only a few terms in the expan- 
sion, then one biases the inversion to give, for example, 
only Gaussian sources. 

We choose to represent the source function in a Basis 
spline (a.k.a. b-spline) basis (^|| as this basis has all of 
the features we require for a good representation. Plots 
of some sample b-splines are shown in Fig. [l] and this 
basis is detailed in Appendix [A]. B-splines are piecewise 
polynomials and are continuous up to the degree of these 
polynomials. The th degree b-spline is the box-spline, 
making our b-spline expansion a natural generalization 
of the approach in Rcfs. |Tl] , |l^ ], Furthermore, in the b- 
splinc basis the concept of the "edge of a bin" in the box- 
spline basis is replaced with the concept of a knot p2[ . 
A knot is simply the place where the polynomials that 
make up the b-spline are patched together. In the "op- 
timized discretization" scheme of Ref. jl^] , the edges of 
the box-splines are varied to minimize the relative error of 
the source. We may generalize this idea to the b-splincs 
easily by varying the locations of the knots. We w ill give 
examples of choosing the knots in subsection [II C and we 
will explain in detail how to choose the "optimal knots" 
in Appendix |a[ 
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C. The reconstruction 

Once we have converted the inversion problem into 
a matrix inversion by choosing a representation of the 
source, we proceed as in Refs. jyjll^Q and extract 
the source. The details of the the Bayesian approach to 
imaging are discussed in the Appendix [b| and we sum- 
marize the main results here. To obtain the coefficients 
of the source, we seek the source that minimizes the x 2 : 



(K • S - lt) T • (A 2 ^)- 1 ■ (K ■ S — TV), 



115) 



where A 2 7?. is the covariance matrix of the correlation 
data. The source that does this is: 

S = A 2 S ■ K T ■ (tfK)- 1 -K (16) 

The covariance matrix of this source is: 



A 2 ^ = {K T • (A 2 ^)- 1 • KY 



(17) 



One should note the dependence on the experimental un- 
certainty A 2 K in Eqs. @ and In Eq. (|l|), the 
points with the largest error contribute to the source de- 
termination the least. Also, in Eq. ([l7]) the points who 
are most effected by the points with the large error also 
have the largest uncertainty. 

In order to stabilize the inversion, we can take advan- 
tage of prior information in the form of equality con- 
straints [poj| . An equality constraint is a condition on 
the vector of source coefficients that has the generic form 
C ■ S = c. One example of such a constraint is that the 
source has slope at the origin. Such a situation arises 
if the normalized particle emission rates, D (c.f Eq. (||)), 
have a maximum. In this case we write 



S'(r 



N M 



0) = 0. 



(18) 



Thus, this case corresponds to Ci = B'^r — > 0) and c = 0. 
We can implement this type of constraint by adding a 
penalty term to the % 2 : 



X 



A(C ■ S - c) 5 



(19) 



Here A is a trade-off parameter and we may vary it in 
order to emphasize stability in the inversion (by making 
A huge) or to emphasize goodness-of-fit (by setting A to 
zero) . Such an ability to trade-off stability for goodness- 
of-fit is discussed in Numerical Recipes Q in detail. 
With this modification of the x 2 , the imaged source is 



S = A 2 S* • (K T • (A 2 ^)- 1 • K + XC T • c) 

and the covariance matrix of source now is 

A 2 S = (K T ■ (A 2 ^)" 1 ■ K + XC T ■ C)' 1 



(20) 



(21) 



There is another way in which prior information enters 
into the inversion - in the representation we choose for 
the source. By using, say, Ni, = 2 in a b-spline expansion 
of the source, we are really assuming that our source and 
its first and second derivatives are all continuous. 

The reader should note that, when we image, we are 
really finding a probability density for the source given 
the correlation data rather than the source itself. The 
set of source coefficients and the covariance matrix of the 
source characterize the height and width of this probabil- 
ity distribution. In the end, we use the source coefficients 
as an estimator of the true source. 



III. TESTS OF THE IMAGING 

We now explore the imaging in the b-spline basis by 
inverting some simple model correlations. We consider 
two model sources, a Gaussian source: 



S(r) 



A 



■ exp (- 



and a source with a dipole form-factor like shape: 



S(r) = A- 



Rd 



(r 2 



4i?i,) 2 ' 



(22) 



(23) 



This second source has a roughly Gaussian peak and an 
extended tail that one could imagine corresponds to long- 
time emission of particles. We chose this source to facil- 
itate comparison to the experimental results in the next 
section. We pick Rq = 4.5 fm, Rd = 3.5 fm and A = 1. 
To generate the correlations, wc convolute the source 
with the proton kernel in Eq. (^|) and bin the correlation 
in 6 McV/c sized g-bins. To simulate realistic data, we 
take the error bars from the real data of Rcf. J27J and add 
statistical scatter consistent with these error bars. The 
data in [^7f are plotted in the same size momentum bins 
as in our test, this data is fully corrected for various ex- 
perimental effects, and these effects are properly reflected 
in their estimates for the experimental uncertainty. The 
resulting correlations are shown in Fig. |^. In all of our 
tests, we confine ourselves to proton correlations because 
the proton FSI arc more important than meson FSI and 
place a more demanding test on the imaging. 

Our first test is to examine how the quality of the 
source reconstruction depends on the number of coeffi- 
cients in the source expansion. In this test we will use 
only box-splines. In the second test, we will use higher 
degree b-splines in the reconstruction. In this test, we 
will use a fixed separation between the knots (equivalent 
to using equal width box-splines). In the third test, we 
will demonstrate the use of the "optimal knots" in anal- 
ogy to the "optimized discretization" method of Rcf. [jl2| . 
In the final test, we will demonstrate the practical use of 
equality constraints. 
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A. Number of coefficients in the source 

Before we begin our imaging tests, we must decide on 
the size of our imaging region and set the number of coef- 
ficients that we wish to reconstruct. Naively, to set r max 
we might use the Fourier estimates from Eqs. ( |Io| ) giving 
Tmax = 103 fm. Experience has shown us that the source 
is usually lost in statistical noise and is consistent with 
zero to within the errors long before this r max . So, we 
set r m ax = 45 fm, roughly half of what the naive Fourier 
estimate suggests. If we find this is too conservative, we 
may increase it later. 

To set the number of coefficients, we could use the 
naive Fourier estimate again. Doing so, we find that 
Ar = 4.1 fm and that we should use 11 coefficients. On 
the other hand, Eq. ( |l3| ) suggests that the imaging prob- 
lem is really a problem of simultaneously solving a set 
of linear equations. Given this, we look at the data and 
see that there are roughly 15-16 points which are differ- 
ent from one and hence contain useful information. This 
suggests that we try using something like 14 coefficients 
in the source. 

This raises the question of the amount of information 
in a data set. If a correlation is Gaussian, than one could 
argue that it contains only two pieces of information: the 
height and width of the Gaussian. On the other hand if 
one bins the data, one could argue that there are really 
N pieces of information corresponding to the number of 
bins where there is an apparent signal. We adopt the 
second viewpoint, but comment that the "amount of in- 
formation" in a data set is an imprecise concept. 

Since it is not clear which number of coefficients we 
should use, we will try both. Additionally, we will image 
using 7 box-splines so that we may compare the results 
with the higher-degree results of the next section. In 
Figs. |^ and || we plot the inversions of the proton corre- 
lations in Fig. ^. 

First we look at the Gaussian source images in Fig. [|. 
In all three panels, the inversions are reasonable represen- 
tations of the true source. Only by looking at the x 2 is it 
clear which image is the "best:" for panel (a), x 2 = 122, 
for panel (b), x 2 = 91 and for panel (c), x 2 — 76. Since 
there are only 83 points in the proton correlations, the 
inversion with 14 coefficients is "too good" and 11 coef- 
ficients seems to be the best choice. Before moving on to 
discuss the dipole sources, we mention that the fluctua- 
tions in the imaged sources are not independent. If they 
were, then we would expect that roughly 1/3 of the bins 
would differ from the true source by at least a standard 
deviation. 

Now we turn to the dipole sources. Looking at Fig. ||, 
none of the images are ideal - all three have large fluctu- 
ations that imply that there is a zero somewhere around 
10 fm. In these three plots, the x 2 is not much of a guide 
either. The x 2 's are 104, 90, and 77 respectively. Finally 



we comment that we cannot tell which images in Figs, y 
and |^ correspond to Gaussian sources or dipole sources. 

We have seen that increasing the number of coefficients 
in the reconstruction helps to reproduce the source, how- 
ever there is a practical limit to how many coefficients we 
may add. As the number of source coefficients increases, 
they become less constrained by the data. At some point 
there are more source coefficients than can be constrained 
by the data and then the extra coefficients only serve to 
reproduce the high frequency statistical fluctuations in 
the correlation data. At this point, the imaged source 
tends to oscillate about the true source as we have over- 
resolved the source. In general, we can never tell when 
we are over-resolving the source as we do not know the 
true source. 

B. Using the basis spline representation 

We expect the source function to be continuous as it 
is the convolution of two emission rates, yet we repre- 
sent it by discontinuous box-splines. Thus, our first im- 
provement over a simple box-spline representation of the 
source is to use higher degree b-splines. 

We re-image the proton correlations in Fig. [2| In 
Figs. | and §, we show the images obtained using the 
b-spline expansion with Nm = 7 coefficients and either 
1 st , 2 nd or 3 rd degree b-splines. Our choice of Nm = 7 
is somewhat arbitrary as we do not know how to extend 
our Fourier transform based estimates to our b-spline ba- 
sis. However, the fact that the imaging works reasonably 
well points to the robustness of the method, despite the 
possibly sub-optimal choice of Nm- In all plots, Nf, + 1 
knots are placed at the end points of the imaging region 
(i.e. at r = fm and r = 45 fm) and the rest of the knots 
are equally spaced between the end points. 

In Fig. |5J, the images are fairly accurate reconstructions 
of the source over two orders of magnitude, but the iVj = 
3 reconstruction is marginally better, due to the higher 
degree of continuity. In all cases, the inversions are better 
than any of the box-spline results. The corresponding 
X 2 's are 99 for degree 1 splines, 94 for degree 2 splines 
and 93 for degree 3 splines. In these plots, the region 
past r = 17 fm is lost in the noise from the correlation. 
We notice that all the plots exhibit the same kind of 
fluctuations seen in the N\, = images in the last section, 
however they are less noticeable because the b-splines are 
so de-localized. Finally, we mention that the unphysical 
rise out past 40 fm is most likely a result of aliasing. It 
is more obvious in these plots because the last b-spline 
has a cusp at the edge of the image. 

We also image the dipole source in Fig. |^(b). Using 
the same settings as for the Gaussian source, we are 
able to reproduce the more complicated behavior of the 
dipole shaped source over two decades in source inten- 
sity. More importantly, upon comparing Figs. Hand pi we 
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can clearly tell the difference between Gaussian and non- 
Gaussian source shapes on the logarithmic scale. This is 
something we could not do with the box-spline represen- 
tation of the source. The x 2 f° r these reconstructions are 
95, 94 and 93 respectively. We comment that the cusp 
in the very first b-splinc helps us represent the relatively 
sharp peak in the dipole source. 

In all of the images shown so far, we see an unphysical 
rise in the source in the far right of the images. This rise 
is most likely a result of aliasing. Chapter 12 of Numeri- 
cal Recipes |2(j] has a detailed discussion of this problem. 
Aliasing is a phenomenon that often occurs when approx- 
imating a Fourier transform over an infinite interval with 
a finite Fourier transform over a finite interval. Consider 
a function f(r) and its Fourier cosine transform, F(q): 

poo 

F(q) = / dr f(r) cos (qr) 

J fr nam (24) 

~ / dr f (r) cos (qr) . 
Jo 

In the first line of this equation, low frequency struc- 
ture in F(q) corresponds to large distance structure in 
/(r), which is neglected in the second line of this equa- 
tion. Now, imagine beginning with F(q) and attempting 
to infer f(r) using a finite r max . What often happens is 
that, whatever strength /(?*) should have out past r max 
gets folded into the region r < r max . In our inversion 
problem, the integral in Eq. (|l|) behaves like a Fourier 
transform. Since statistical fluctuations in the data are 
artificial high-frequency structure, we should not be sur- 
prised to see features reminiscent of aliasing when we 
image. Based on our experience, adjusting r max or con- 
straining the source at r max can help cure this problem. 
However, the rise at the largest r is usually preceded by a 
region of the image that is consistent with zero so we can 
easily identify the usable part of the image and ignore 
any artifact due to aliasing. 

C. Choosing the knots 

For our next refinement, we examine how choosing the 
knots affects the inversion. Were the problem of imag- 
ing as simple as inverting a Fourier transform, the op- 
timal bins in r would be evenly spaced and given by 
Eq. (|l0|). However, the kernels we are interested in arc 
often distorted by the Coulomb repulsion of the pairs as 
well as other FSI. Furthermore, some regions of the data 
have large errors and it would be useful if we could com- 
bine those bins somehow. Taken together, we must ask 
whether keeping equally spaced bins in the source is op- 
timal. In Ref. 0, we found that we could improve the 
imaging dramatically by choosing the size of the bins in 
r to minimize the error in the source relative to a test 
source. This technique can be generalized to the b-splinc 



basis by simply varying the knots. To choose the "opti- 
mal knots" we proceed as mentioned earlier and detailed 
in Appendix |A[ 

For the Gaussian source in Fig. [t], wc show the in- 
versions using 3 rd degree b-splines using 7 coefficients 
for both fixed knots (as in Fig. ||(c)) and the "optimal 
knots." Several things are apparent in this figure. First, 
the fixed width knot reconstruction is markedly better 
than the lower-degree b-spline reconstructions in the pre- 
vious section, simply due to the higher degree of conti- 
nuity. The x 2 °f this reconstruction is 93. Using the 
"optimal knot" reconstruction, the source is everywhere 
consistent with the true source except at the lowest r 
(< 5 fm) where the source drops nearly an order of mag- 
nitude. This drop is unphysical for a source which is the 
convolution of two single particle sources, each with one 
or more maxima. This drop is due to the close pack- 
ing of the knots at the lowest r and can be remedied by 
lowering the number of coefficients in the reconstruction, 
increasing the size of the fit region, or by using equality 
constraints (as we will show in the next section). The x 2 
of this reconstruction is 90. 

In Fig. ||, we show the similar set of inversions for the 
dipole shaped source. Both inversions seem to do a rea- 
sonable job of representing the source, except at the low- 
est r where the cusp of the first b-splinc is a bit higher 
than the true source. The "optimal knot" reconstruction 
is marginally better that the equally spaced knot recon- 
struction as it has a % 2 of 92 compared to a x 2 of 93 for 
the equally spaced knot reconstruction. 

Given the inconsistent performance of the "optimal 
knot" reconstruction, we ask ourselves why this refine- 
ment does not always help. To find the optimal knots, 
we move the knots to minimize the error on the source 
relative to a test source (which is the same for all of the 
inversions). The error on the source depends only on the 
kernel and the error in the data so the "optimal knots" 
do not know anything about the true source. If the true 
source has interesting structure someplace where we are 
not sensitive to it, the "optimal knots" will be widely 
spaced there and we will not have the resolution to see 
the structure. Conversely, the "optimal knots" may end 
up giving very high resolution exactly where we do not 
need it, witness Fig. 0(b). On the other hand, the "opti- 
mal knots" can give resolution exactly where we need it, 
as in Fig. ||(b). 

D. Equality Constraints 

Now we come to the final refinement, the use of equal- 
ity constraints. As we have mentioned before, a con- 
straint is a piece of prior information such as knowing 
that the first derivative at r = should vanish for a dif- 
fercntiable angle-averaged source, S'(r — > 0) = 0. Using 
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constraints amounts to adding information, so we imag- 
ine that we will be able to use more coefficients in the 
reconstructions. This we will see illustrated below. A list 
of constraints we could use are shown in Table ffl. 

In Fig. ||, we show inversions using the S'(r — > 0) = 
constraint for the Gaussian source. We used the "optimal 
knots," 3 rd degree b-splines and 7 coefficients (in panel 
(a)) and 9 coefficients (in panel (b)) in these inversions. 
We see that we have solved the pathological behavior of 
the imaged source at low r and the agreement with the 
true source appears good. Upon examining the x 2 (107 
for panel (a) and 93 for panel (b)) we see that the 7 
coefficient source in a lot worse than it appears as it is 
routinely higher than the true source in the region from 
10-20 fm. In Fig. [l(], we show the results of using the 
same constraint to image the dipole-shaped source. Here 
we see that, for 7 coefficients (panel(a)), the quality of the 
image has gone down considerably: we no longer match 
the height of the peak and we cannot resolve any of the 
tail. The x 2 for this inversion is an comparatively large 
108. For 9 coefficients (panel(b)), the situation is much 
better. We now get the peak and can resolve the tail. 
The x 2 nere is 89. 

We see that this one constraint gave us the ability to 
add another two points in the reconstructions without 
over-resolving the source. At a practical level, the first 
few b-spline coefficients must be adjusted together in or- 
der to satisfy the constraint, in effect leaving fewer co- 
efficients to fit the data. Thus, we must add more co- 
efficients if we want to simultaneously fit the data and 
satisfy the constraint. This observation fact leads us to 
posit a rule of thumb: "amount of information in the 
data" + "number constraints" > "number of coefficients 
in expansion." Additionally, one should pick the number 
of coefficients somewhere near what one would estimate 
based on the Fourier estimates discussed earlier. 

Finally, by introducing all three refinements of the 
imaging (b-splines, "optimal knots," and equality con- 
straints) we are able to reproduce the height of the source 
at r = quite accurately. The value of the source at r = 
is essential for extracting the space-averaged phase-space 
density pi| , P8|| . 

IV. ANALYSIS OF NA44 CORRELATIONS 

Since we can reliably image a source from correlations 
using the Bayesian approach to imaging in a b-spline rep- 
resentation, we turn to the analysis of NA44 correlations. 
In a series of papers HU, NA44 detail their measure- 
ments of angle-averaged pion, kaon, and proton correla- 
tions from p+Pb collisions at 450 AGeV/c and central 
S+Pb collisions at 200 AGcV/c. In two of the earlier 
papers [f^H, they claim to have detected non-Gaussian 
kaon and pion correlations. To bolster their claim, they 
fit the Coulomb corrected correlations to Gaussians and 



to exponentials. In particular, they fit to the following 
functional forms: 

1l{Q mv )=\cxv{-Q 2 nv R 2 G ) (25) 

implying the Gaussian source of Eq. (|22| ) and 

K(Qinv) = Acxp (-2Q mv R D ) (26) 

implying the source with a dipolc-like shape of Eq. (|2^). 
Here Q m „ = 2qi nv = y— (pi — P2) 2 , the relative momen- 
tum variable traditionally used in the analysis of meson 
correlations. The NA44 correlations that we image are 
collected in Fig. |ll|. 

In this section, we first image the NA44 correlations. 
Second, we compare the images to the results of some of 
NA44's fits. Next, we discuss NA44's RQMD simulations 
of the S-Pb reaction and the implications for the source 
function images. Finally, we discuss the sources from the 
NA44 p-Pb data. 

A. Imaging Analysis 

The results of the imaging analysis are presented in 
Fig. [l2|. As a cross-check, in Fig. [ll] we plot the corre- 
lations corresponding to the inverted sources along with 
the original data. In these inversions, we used either the 
noninteracting meson kernel in (0) (for the Coulomb cor- 
rected pion and kaon correlations) or the proton kernel 
in (^J). Due to the differences in kernels, binning and 
quality of the various data sets, each image had to be 
hand tuned separately. Since we do not know the true 
sources that correspond to the data, we used a set of 
three criteria to decide when we have a good source: 

1. Is the image stable - i.e. when we tweak a pa- 
rameter (e.g. the number of bins, r max , number of 
constraints, etc.) does it change much? 

2. Does the imaged source give a correlation consis- 
tent with the original? 

3. Is the x 2 as small as we can make it? 

In all cases, we used third degree b-splines. The param- 
eters of the inversions are collected in Table O. Only the 
p-Pb pp source was imaged without the r = smooth- 
ness constraint. We did not use this constraint because 
the knot density is too low at low r. Turning on the 
constraint widens the peak artificially as the next few 
b-splines have to be tuned to get the zero slope at the 
origin, dramatically raising the x 2 - 

When looking at the images, several things are appar- 
ent. First, each of the sources from the p-Pb reactions 
are roughly a factor of two narrower than the correspond- 
ing sources from the S-Pb reaction. This is most likely 
a result of the different system sizes. Second, a com- 
parison of the sources from the same reaction reveals 
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that the pion sources are wider than the kaon sources 
and the kaon sources are wider than the proton sources. 
Next, it is apparent that all six of the sources have main 
peaks that appear Gaussian. However, both the pion and 
kaon sources in the S-Pb reaction have significant non- 
Gaussian tails. These tails are most likely not due to 
aliasing as r max in both plots is at 35 fm, but, in order 
to show all six sources on the same scale we truncated 
the plots at 20 fm. Unfortunately, this means that we 
cannot display that the kaon and pion sources are con- 
sistent with zero in the region from 25-30 fm nor can one 
see the rise that is obviously due to aliasing past the re- 
gion 30-35 fm in the kaon source. No aliasing is apparent 
in the pion source or any of the proton sources. In the 
pion and kaon sources from the p-Pb reaction, aliasing is 
apparent on the far right side of the plots. 



B. NA44 Fits 

Following the imaging analysis of the NA44 correla- 
tions, it is useful to compare those results to the various 
fits performed by NA44. We have in mind two sets of 
fits: the one-dimensional fits to Gaussians and exponen- 
tials in Rcf. [QJ|] and the three-dimensional Gaussian fits 
in Rcf. HJ. 

We first compare the imaged sources to the results of 
NA44's one-dimensional fits. In addition to the imaged 
sources, in Fig. [l| we show the fits as solid curves (for 
the Gaussian fits) or dashed curves (for the exponential 
fits). NA44's fit parameters are summarized in Tables III 
and IV. In the S-Pb sources in Fig. O, both the kaon 



and pion sources seem to be consistent with both the 
Gaussian and the dipole-shapcd curves in the range from 
4 fm out to about 12 fm. At the lowest r the pion image 
seems to split the difference between the two fits and 
the kaon image is below both fits. Both sources exhibit 
long non-Gaussian tails. In the pion source, this tail is 
higher than the dipole fit and in the kaon source, the tail 
is consistent with the dipole fit. For the p-Pb sources 
in Fig. both the kaon and pion sources appear very 
nearly Gaussian. At the lowest r, the kaon source is a 
bit below the Gaussian fit. 

One might ask if the tails in the sources are sim- 
ply due to the non-spherical geometry of the full three- 
dimensional source. To answer this question, we com- 
pare the imaged sources to sources corresponding to 
NA44's three-dimensional fits in Ref. |3(J. These cor- 
relations were acquired as follows: NA44 first measured 
like-charged pion and kaon correlations in the Longitudi- 
nal Co-Moving System (LCMS), Coulomb corrected their 
data, and then fit their correlations to a Gaussian form. 
The form they fit to is: 



C(Q) = 1 + Aexp 



Q L R L 



The results of these fits are shown in Table [v| for 7r + 's 
and K +, s. Now, as in the one dimensional case, a 
three-dimensional Gaussian correlation corresponds to a 
three-dimensional Gaussian source. For the correlation 
in Eq. (p7[), the corresponding source is 

S(r) = 



{2^fR Ts R To R L 



x cxp 



1 



'Ts 



4\R Ts 



'To 

R 2 

n To 



it''' ^ 



Looking at the various transverse fit radii obtained by 
NA44, it seems safe to set Rt s = Rto = Rt m this 
source. To convert this three-dimensional source to a 
one-dimensional source, we only need to angle-average 
it, S(r) = ^J dn r S(r), to obtain 

^2 



S(r) 



A 



16nR T R L 

2R eff 



cxp 



x 



r 

2R eff 



erf 



crfi 



4R T 



2R 



eff 

r 



for R T > Rl 



2R, 



(29) 



for R T < Rl 



(27) 



where R e ff = Vy \Rl 2 ~ Rt' 2 \ an< ^ cr ^( x ) = ierf(— ix). 
With this result, we are now able to compare the imaged 
sources to the angle-averaged three-dimensional fits from 
NA44. This is shown in Fig. 

For these fits, we used A and Rl as given Table and 
Rt was chosen to be the average of Rt s and Rt ■ Exam- 
ining Fig. [ll|, it is clear that despite the fact that NA44's 
three-dimensional sources correspond to non-Gaussian 
angle-averaged sources, they are not able to account for 
the extended tails we find in the images. There is a 
caveat with this conclusion: their fits were performed 
in the LCMS and our source is imaged in the pair CM 
frame. Since the two frames coincide only in the limit 
that pt — > 0, we computed the pion source using the 
low-pT data set only. Only one data set is available for 
the kaons. Nevertheless, it seems unlikely that the pion 
or kaon fit parameters would change dramatically if the 
fits were performed in the pair CM. 



C. Discussion of S-Pb Data 

In addition to the various fits to the correlations, mem- 
bers of the NA44 collaboration performed RQMD simula- 
tions §-|,|l3],||] of the S-Pb and p-Pb collisions. Rather 
than repeat this work, we will summarize it and explain 
its implications for the source shape. In all but the pp 
case, the simulated RQMD correlations compared favor- 
ably with the data. In the pp case, the simulations over- 
estimate the height of the correlation peak by roughly 
30%. 

Sullivan et al. Jl3[ | explain that the width of the kaon 
correlation (corresponding to the width of the imaged 
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kaon sources) is determined mainly by the size of the 
kaon production region. First, kaons are mainly pro- 
duced directly in the reaction (either from fragmenting 
strings or hadronic reactions) or from the decay of K* 
resonances. Now, the reaction zone is roughly the size of 
the sulfur nucleus (R rms =3.3 fm) and this defines the 
size of the Gaussian core. The K*'s are also produced 
in a region of roughly this size. However, since their life- 
time is roughly r w 4 fm/c, K*'s do not travel far before 
decaying, giving rise to a non-Gaussian halo surrounding 
the Gaussian core. This halo is neither exponential nor 
Gaussian, but rather a convolution of a Gaussian source 
with an exponential []l4| , [l5| . Since the kaons have a much 
larger mean-free path in nuclear matter than cither pions 
or nucleons (at least in RQMD), they do not rescatter as 
the system evolves. This means that their last interaction 
point is very nearly the size of this production zone. 

The pion width in the S-Pb reaction is also described 
well by the RQMD calculations. The initially pro- 
duced pions (produced mainly from hadronic reactions, 
although there is a string component) also have a source 
width set by a combination of the geometrical overlap 
of the colliding nuclei and subsequent dynamics of the 
system. From their production until they freeze-out, 
the pions interact strongly with the system. Because 
this evolution involves longitudinal flow, there arc strong 
position-momentum correlations in the freeze-out posi- 
tions of the pions. The longitudinal size of the region 
where one can find pions with a low relative momen- 
tum is comparable to the transverse size of the system at 
freeze-out giving a Gaussian core somewhat larger than 
the initial system size. In addition to this core, the pions 
also have a large contribution from the decay of various 
resonances, mainly the p, u>, rf, and 77°. Now the the 77's 
are not capable of altering the pion source size or shape 
as the 77's lifetimes are much too long. On the other hand, 
both the p (whose lifetime is f» 2 fm/c) and the ui (whose 
lifetime is ~ 23 fm/c) both contribute to a non-Gaussian 
halo in a manner analogous to the K* 's above. Since the 
p'a lifetime is smaller than the u>'s, it contributes to the 
shorter distance part of the tail and the u> to the longer 
distance part. 

In the case of the protons in the S-Pb reaction, the size 
of the source is mainly set by the size of the region where 
we find protons with low-relative momentum at freeze- 
out. Looking at the plots of RQMD correlations in [p9[ , 
the RQMD simulations are somehow unable to reproduce 
the sizes of these regions. Since the size of this region is 
a direct function of the geometry of the source and the 
space-momentum correlations in the source (these cor- 
relations are caused primarily by flow), RQMD's failure 
here is somewhat of a mystery. Unfortunately, our image 
is not detailed enough to give much of a clue where the 
discrepancy arises. A higher resolution measurement of 
this correlation would help immensely 



D. Discussion of p-Pb Data 

Unlike the S-Pb reaction, RQMD is able to describe all 
of NA44's measured p-Pb correlations quite well [7|,p|,p9|. 
Because the p-Pb system so is small, the reaction is most 
likely dominated by the formation of a few color strings 
and/or ropes. RQMD uses the Lund string model to 
model string formation and fragmentation |3l| ] , so to un- 
derstand what is happening in the sources it is useful to 
understand some of the features of this type of model. 
In the Lund model, momentum and spatial rapidities are 
tightly correlated. We will show that this correlation 
leads to an approximate scaling of the K and p source 
radii with mass: 

fnxRx ~ rripRp. (30) 

As we will point out, this scaling is born out by the im- 
ages. We mention that the pion source should not follow 
this scaling because a large fraction of the pions result 
from resonance decays which distort the shape of the pion 
source. 

In a Lund-type string model ]32] |, a meson is viewed 
as a qq pair attached by a string and this pair oscillates 
back and forth in the linear confining potential of the 
sting. This type of model can also describe baryons if one 
replaces one of the quarks at the end of the string with 
a di-quark. In this discussion, we ignore the transverse 
extent of the string and imagine the string lives in a 1+1 
dimensional space. As the qq pair oscillates, in one period 
it sweeps out an invariant area given by 

2 

m 

Ax + Ax_ = — (31) 

in light-cone coordinates. Here m is the mass of the 
hadron and a is the string tension. A hadron breaking 
off from a string is pictured in Fig. |^| 

To assess any correlations between the mass and the 
production location of a hadron, we examine the distri- 
bution of break-up points of the string that produces the 
hadron. In Fig. [u], points 1 and 2 are the locations where 
a qq pair separates from a string. Assuming a constant 
break-up probability per unit time per unit length, P, 
the distribution of break-up points is given by 

dV oc exp {-Px W xf)6{x W - a£ ) )0O B + ) " x + } )- 

(32) 

The exponential in this expression gives the probability 
that the string does not break up before the q and q are 
made at points 1 and 2 and the theta functions ensure 
the proper ordering of the coordinates. We define the 
creation point of the hadron of interest as the mean of 
the break-off points of the qq pair: 

X± = \{x ( i ) +xf). (33) 



9 



We also write 

Ax± = ±{x y ±> - 
In terms of these, we may write Eq. ([52]) as 
dP oc 0(Ax_)0(Aa;+) 

1 



(34) 



x exp I -P yX+X- + -Ax + Ax_ 
+ ^(Ax^X+ + Ax+X_) 



(35) 



For a hadron at mid-rapidity, from Eq. ( |3l"| ) we have 
Ax + = Aa;_ = -j=^- Writing the hadron position back 
in Cartesian coordinates, wc find 



dP oc exp 



T -X z 



T )\9{T-\X\) 



(36) 



This last theta function makes causality explicit. 

Since dP expresses the probability of creating a hadron 
at position X and at time T, it is proportional to the 
emission rate in Eq. (0). Thus, we can imagine doing the 
convolution in Eq. (@) to obtain the source size. Given 
that the emission rates in this case are not Gaussian, we 
should not expect the source function itself to be Gaus- 
sian. Nevertheless, we can still estimate the width of the 
source function from Eq. (|36"|). We take the source width 
to be the distance at which the magnitude of the source 
function drops by 1/e. From Eq. (p6|), wc sec that the 
probability of creating a hadron drops by 1 /e by roughly 



T ~ X 



2a 
mP' 



(37) 



This implies that the source function itself will have a 
width which is correlated with the mass of the hadron 
and given by 



I? 



2\/2cr 
mP ' 



(38) 



Here, the factor of \f~2~ arises from the convolution in 
Eq. (H). If wc make the reasonable assumption that the 
string tension and break-up probability are both uni- 
versal, than we have the scaling relation in Eq. (30). 
Looking at the kaon source in Fig. |l|(d), the kaon 
source 1/e- width is Rk ~ 6 fm. Following the scal- 
ing, the proton source should drop by 1/e by roughly 
R p « ^^Rk = 3 fm. This is roughly born out in the 
image in Fig |l^(f). With (|38|), we can go further than 
than just checking this scaling and try to compute the 
source width directly. Taking typical values for the string 
break-up probability, P = 1 fm 2 , and for the string ten- 
sion, cr=l GeV/fm, we find that the proton source 
radius should be about R=2.9 fm and the kaon source 
radius should be about R=5.7 fm, again in rough agree- 
ment with the images. It would be very interesting to 



see if this scaling persists in sources imaged from other 
like-pair correlation in p-A collisions such as pp, KqKq 
or A A. 

In this consideration, we have neglected several things. 
First, although we have ignored the transverse degrees of 
freedom of the strings, the longitudinal length of the pro- 
duction region is much larger than the transverse extent 
so it is the longitudinal extent that determines the angle- 
averaged source radius. Second, we have ignored the 
finite mass of the quarks. Adding finite quark masses 
would change the trajectories of the qq pair and the 
shape of the area swept out by the pair, but it should 
not change our conclusions appreciably. Third, we have 
neglected any consideration of other strings in the p-Pb 
interaction zone. Next, if there are more than one string 
produced in the reaction, the possibility for string fu- 
sion exists and it is neglected here. Finally, we have ne- 
glected the Bjorken-like position-momentum correlations 
along the string length. These correlations will narrow 
the source |]l6| and possibly account for the minor differ- 
ence between the predicted and measured proton source 
widths above. 



V. CONCLUSION 

In this paper, we have investigated the possibility of de- 
tecting non-Gaussian sources in heavy-ion collisions. In 
simple, but realistic, model calculations we have demon- 
strated that it is possible to distinguish between Gaussian 
and non-Gaussian source shapes using an improved imag- 
ing method and high resolution data. Imaging not only 
has achieved results comparable to Gaussian fits, but has 
now uncovered deviations from Gaussian behavior. 

This improved imaging method has several features 
that make it superior to the previous methods in pl| , ^2l . 
First, this method uses Basis splines to represent the 
source, giving a continuous representation of the source. 
The resolution of the images is controlled by the place- 
ment of the knots (the points where the polynomials that 
make up the spline basis are matched). Since the knots 
are analogous to the edges of the bins of a box spline rep- 
resentation, we may do as in [fl2f and use the "optimal 
knots" to further improve the image. In addition to these 
improvements, equality constraints are now implemented 
in a simple manner. As in the previous imaging method, 
imaging is still a least-square problem in which the co- 
efficients of the Basis spline representation are chosen to 
minimize the \ 2 of the data. In other words, the imaging 
gives the source that has the highest probability of rep- 
resenting the correlation data. Finally, we mention that 
the amount of information available in the data limits the 
amount of information we may extract in form of image 
and constraints "increase" the amount of information. 

Using this improved imaging method, wc analyzed the 
proton, kaon, and pion correlations from S-Pb and p- 
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Pb reactions measured by NA44. We find evidence for 
non-Gaussian halos in the source functions from the like- 
meson correlations in the S-Pb reaction. These halos are 
likely to be due to resonances decaying and producing 
the mesons. RQMD model simulations of the correlations 
reproduce the experimental data quite well except in the 
case of the proton correlations from the S-Pb reaction. 
Unfortunately the source image does not shed much light 
on the discrepancy. In the case of the p-Pb reaction, the 
imaged sources suggest a scaling of the source widths 
that we should expect on the basis of Lund-type string 
models. This scaling should be tested by examining other 
like-pair correlations (e.g. pp and AA) where these other 
pairs do not suffer from large resonance contributions. 
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APPENDIX A: BASIS SPLINE 
REPRESENTATION OF SOURCES 



In this paper and previously in 1 12 1, we write the source 
expanded in some basis: 



N M 



(Al) 



In |Q this basis was the box-spline basis. In the box- 
spline basis, the widths of the individual bins could be 
varied to increase or decrease the resolution of the kernel 
to minimize the relative error of the source. Unfortu- 
nately, in this representation the source is not continu- 
ous. 

In this work, we expand our sources in a more general 
basis: that of Basis splines, a.k.a. b-splines |22]. B- 
splines are piece-wise continuous polynomials which can 
be made arbitrarily smooth by changing the degree of 
the splines; the th degree b-splines are actually box- 
splines used before. The b-splines are characterized by 
a set of knots which mark the points where the various 
polynomials that make up the b-spline are matched. In a 
sense, these knots generalize the "edges of the bins" of the 



box-splincs. For this reason, we may vary the locations 
of the knots to minimize the relative error of the source, 
generalizing the method in |l2[ . 



1. The basis splines 

Now we define the b-splines. A N£ h degree b-spline is 
characterized by a set of knots, {tk} with t\ < t-2 < ... < 
tN knots - The number of knots must be chosen so that 
NkZL > N M + N b + 1 for Nm b-splines. For N b = 0, 
the b-splincs arc box splines, i.e. 



Bo,i(r) = Xi(r) 



1 if tj < r < t i+1 
otherwise. 



(A2) 



Note, if ti = ti + i then Bq^(t) — 0. The rest of the b- 
splines may be constructed from this first one through a 
set of recurrence relations 



B Nb+ i,i{r) = w Nb+1A (r)B Nbti (r) 

+ (1 - w Nb+1 , l+1 {r))B Nb , l+1 (r) 



(A3) 



where the weight factor is 



w Nb +i,i( r ) = { ti+N b 




if ti ^ t i+ N b 
otherwise. (A4) 



With these definitions, we may write a b-spline of any 
degree back in terms of the box-splincs: 



i+N b -l 

B N b A r ) = Yl h N„-3 X A r ) 



(A5) 



where 6jv bJ is a polynomial in r of degree N b which we 
will not write explicitly. We plot sample b-splines of dif- 
ferent degrees in Fig. [y. 

There are three other properties of the b-splines that 
are of note. First, b-splincs arc normalized so that 



drBi{r) 



i+Nb+l ~ h 

N b + 1 



(A6) 



Second, the i th b-spline is zero outside of the region ti < 
T < ij+iVb+i. Third, the b-splincs arc not orthogonal but 
expressions for their inner product exist ]22| . 

From the definitions and from the figure, it is not clear 
how to pick the knots. The knots do not have to be 
equally spaced and, in many situations, it is best not to 
space them equally. In fact, one can even pile up to N b +l 
knots on the same point! One can do this because 

N b + 1 = number knots at a point + number continuity 
conditions at that point. 
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If a b-splinc has N^ + l knots at a point then that b-splinc 
is discontinuous there. Away from these multiple knots, 
the b-spline is still continuous up to derivatives of degree 
N . In the main text, we remove all assumptions about 
the continuity of the source at the boundaries by keeping 
Nt, + 1 knots at the boundaries. We then optionally re- 
insert continuity conditions using equality constraints. 



2. Optimizing resolution 

Now we come to the point of choosing the best knots 
for the inversion, generalizing the "optimized discretiza- 
tion" scheme of Rcf. Jl2l. First, the model covariancc 



matrix (cf. Eq. (B7)) depends on the kernel of the inver- 
sion, the error on the data and whatever scheme we use 
to represent the source, but not on the data or the model 
itself. For a given kernel and set of data errors, we are 
free to change our representation of the source in order to 
minimize the error of the source. In particular, we may 
vary the location of the knots (at least not the knots fixed 
at the endpoints of the imaging region) to minimize the 
error of the source coefficients, ASj = -J A 2 Sjj, relative 
to some dummy source: 



N knots -N b -l 

E 

j=N b +2 



AS,- 



S 



dummy 



(A7) 



The coefficients S^ ummy are the coefficients of the ex- 
pansion of a dummy source in b-splines. In this mini- 
mization, the first and last + 1 knots are held fixed 
and the positions of all of the other knots are varied. 
The dummy source itself is chosen to be big roughly 
where one expects the source to be big and small where 
one expects the source to be small. Since the detailed 
shape of the dummy source should not be important, 
in this paper we chose an exponential dummy source 
with radius R dumm v = 3.5 fm given by S dummy (r) oc 
exp(-r/i? d "" ims/ ). 



APPENDIX B: BAYESIAN APPROACH TO 
IMAGING 

In this appendix, we will explain the technical details 
of the Bayesian Approach to imaging and extend the 
approach of Ref. Jl2| by implementing constraints in a 
more consistent manner. In the previous works the con- 
straints are implemented by Monte-Carlo sampling the 
experimental errors, leading to statistical fluctuations in 
the extracted source. In that approach, no distinction 
is made between equality and inequality constraints. By 
equality constraints, we mean constraints of the form 



dr/(r)5p(r) = constant, 



(Bl) 



and by inequality constraints, we mean constraints of the 
form 

J dr/(r)S* P (r) > constant. (B2) 

Both of these types of constraints are forms of prior in- 
formation, meaning information we have in hand before 
we began imaging. In this appendix, we will explain how 
to use prior information in an inversion. In particular, we 
will discuss two methods for implementing equality con- 
straints directly into the inversion and we will discuss a 
better way to implement inequality constraints. We con- 
clude this appendix with a brief discussion of inequality 
constraints. 



1. General Theory 

Suppose we have an Nu dimensional vector of observed 
data d ofcs with covariancc matrix A 2 d that we wish to 
represent by some iVjvf dimensional model vector m. In 
the case of source imaging, our model vector is the vec- 
tor of coefficients of the function expansion of the source 
and the data is the raw correlation function. We assume 
the data has a diagonal covariancc matrix. In an ideal 
measurement of the data, there would be no experimental 
uncertainty and the data and the model would be related 
through some linear equation: 



d = K ■ m. 



(B3) 



Here, K is the kernel of an integral equation. 

In a real imaging problem, the observed data has er- 
rors and statistical scatter. To make progress, we adopt 
the so-called Bayesian Approach to imaging where we 
seek the probability density, er(m), for a specific model 
m to represent the data p3p3| . With this density, we 
take the mean of the density as an estimator of the true 
model and the width of the density as an estimate of the 
uncertainty in our model. Neglecting the error in our de- 
termination of the kernel and assuming the uncertainty 
in our measurement of d ofcs is Gaussian, we can write 
Bayes Theorem as follows: 



er(m) oc p(m) exp 



1 



"Kdata 



(B4) 



Here, p(m) encodes all of the prior information we have 
about the model and Xdata 

, T 



Xdata =(K-m 



iobs \ 



{A 2 d)~ 1 (Km-d 



obs > 



(B5) 



Here the superscript T represents a matrix transpose. 
The dimension of the model vector, Nm, and the dimen- 
sion of the data vector, No, need not be equal. Indeed, 
it is better to have many more data points than model 
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parameters so that we may over-constrain the system. 
For the time being, assume we have no prior information 
so we may set p(m) = 1. 

We immediately see that the most probable model vec- 
tor is the one that maximizes the probability and hence 
minimizes the xLta- Following we 
can find both the model vector that minimizes this Xdata 
as well as the covariance matrix of the model: 



where 



(m) 



A 2 m • K T 



(A 2 dy 



obs 



and 



A 2 m 



(K 1 



(A 2 ^- 1 -K)-\ 



(B6) 



(B7) 



Eq. (B6) usually goes by the name of a normal equation. 
The model covariance matrix is independent of d obs and 
depends only on the error of the data and the kernel 
itself. 

We can write the Xdata directly in terms of the model 
covariance matrix and the Xdata minimizing model vec- 
tor §3: 



Xdata = (m- <m)) J 



(A 2 ™)" 1 • (m- (m)). 



(B8) 



In other words, if we have a Nm dimensional model space, 
then the Xdata = 1 hypersurface is an A^M-dimensional 
hyper-ellipsoid with principal axes given by the eigenvec- 
tors of the covariance matrix of the model, A 2 m. These 
principal axes do not need to correspond to the directions 
corresponding to the nij components in the model space. 



2. Equality Constraints 

Now let us discuss the role of prior information in the 
general inversion proble m. We concentrate equality con- 
straints such as in Eq. (Bl). In a matrix form, equality 
constraints are written as 



C ■ m = c. 



(B9) 



Here C is a matrix of constraint equations and c is a con- 
stant vector of constraint values. In the inversion prob- 
lem in the main text, there are a variety of constraints 
we might use and they are listed in Table |. The equality 



constraint in (B9) corresponds to the prior probability 
density of 



p(m) oc S (C ■ m — c) - 



(BIO) 



Equality constraints can be cast into a Gaussian prior 
probability density simply by writing this density as a 
Gaussian with vanishing width: 



-\xl q uai 



(Bll) 



xlqual = (C ■ m - C)' 



(B12) 



Finding the most probable model then corresponds to 
minimizing a modified % 2 : 



~"2 _ 2 , \ -A 
X Xdata ' "Xequal • 



(B13) 



The solution is straightforward and corresponds to the 
most probable model: 



i) = A 2 m • (K T ■ (A 2 **)" 1 • d° bs + \C T 



along with the model covariance matrix: 

A 2 m = (K T ■ (A 2 ^- 1 ■ K + XC T ■ C) 



(B14) 



(B15) 



It is clear that to correctly simulate a delta func- 
tion prior information density, we must choose a large 
A. Looking at Eq. (B13), we must choose a A so that 

xlqual > Xdata ' We n0W estimate the Sizes Of xlqual 

and Xdata- A good fit to the data should have the Xdata 
nearly at the number of degrees of freedom, i.e. Xd„ t„ ~ 
N D - N M - To estimate xl qua u we follow Ec l- (|BlSj) . In 
the main text, a typical source is > 1 x 10~ 6 fm~ 3 and 
the constraint matrix and vectors are typically ~ 1 fm 3 , 
so xlquai ~ Nm x 10~ 12 . Putting this together, we must 
choose 



^ ^ > Xdata/ ' X equal 



(Nn 
\N m 



x 10 1 



(B16) 



For example, for Njj ~ 83 and Nm = 8 we need A ^> 
10 11 . By adjusting strength of A, we can adjust strength 
of various terms, emphasizing stability of inversion (i.e., 
obeying constraints) over representing the data. Thus, A 
functions as a trade-off parameter in the jargon of inverse 
theory. See section 18.4 of Numerical Recipes |pr| for a 
more complete discussion. 

A useful alternative to this scheme (and a way to do 
A — > oo limit exactly) is to use the Householder transfor- 
mation to eliminate the constraints from the unmodified 
normal equations of Eqs. (B6) p^J 



. The trade-off is 
that the Householder transformation may be somewhat 
unforgiving. Due to an unfortunate choice of basis func- 
tions, it may not be possible to satisfy two constraints 
simultaneously even if they can be satisfied simultane- 
ously in the true answer. By keeping A finite, we are 
never trying to satisfy the constraints exactly so can do 
a reasonable job of obeying both constraints. Neverthe- 
less, schemes based on Householder reductions of the con- 
straints are complementary to ones using the Gaussian 
prior probability. 
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3. Inequality Constraints 

Now we ask how to use constraints of the form 

Cm>c. (B17) 

Such constraints are called inequality constraints and 
there are many different ones we could use: the source is 
positive definite, the derivative of the source is bounded 
(to ensure smoothness), or the source satisfies the Fourier 
transform test from Ref . |Tc| ] . Inequality constraints cor- 
respond to prior probability densities of the form 



p(m) oc 9 (C ■ m — c) 



(B18) 



which cannot be rendered into a Gaussian form. 

In Ref. pi , we use a simple Monte-Carlo sampling 
scheme to implement inequality constraints. In this 
scheme, one uses the experimental uncertainty to gener- 
ate an ensemble of correlations, each consistent with the 
original. One then inverts each one to obtain a sample 
source and discards any sources that are not consistent 
with the inequality constraints. One then combines the 
samples that are consistent with the constraints to obtain 
an average source and an estimate of the errors on the 
source. The problem with this scheme is that it pushes 
the sources away from edges of the model space defined 
by the constraints. 

We illustrate this problem with a simple example. Sup- 
pose we have an inversion problem where the goal is to 
determine two points S± and 5*2 under the constraint that 
52 > 0. We sketch one possible outcome of the inversion 
in Fig. [15[ In this picture, we see the best fit value of 
Si and S*2 is consistent with our inequality constraint, 
but the constraint cuts through both the la and 2a con- 
tours. Using the Monte-Carlo sampling scheme discussed 
above, we would actually be finding a false best-fit point 
which is slightly above and to the left of the true best-fit 
value because we throw out samples with S*2 < 0. The 
errors on these points would also be symmetrically placed 
around this point. In fact the correct way to solve the 
problem is just to quote the best-fit values of S\ and 52, 
with asymmetrical errors. 

The way inequality constraints are implemented in 
most commercial inversion packages is through so-called 
"Active Set Methods" f35j . In these methods, one finds 
the best-fit solution as one normally would have if there 
were no inequality constraints. If the best-fit solution 
lies in a region excluded by the inequality constraints, 
then the code finds the edge of the included region (the 
so-called Active Set) and searches along it until the code 
finds the solution that minimizes the x 2 . Such a scheme 
is powerful, but likely beyond what is needed for our 
problem. 
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TABLES 



constraint 


continuous representation 


b-spline representation 


flat at r — 


£(r->0)=0 
or 


^S f S^(r^0) = 

3=1 


normalized to A 


A— I J~ ^,2 Q/ \ \ 

47T / (IT T 0\T) — A 

Jo 


N M poo 
> 47TOj / fir 7 Dj\T ) — A 
3 = 1 ^° 


zero outside of imaged region 




5^5 J B,-(r m « H ,) = 

i=i 


flat at T = Tmax 


a \ ' max } — ^ 

or 


^ SjBj[r m ax) = 



TABLE I. Equality constraints on the b-spline representation of spherically symmetric sources. 







Tmax [fm] 


# coeffs 


r — constraint? 


# data pts. 


x 2 


S-Pb 


TT + 7V+ 


35 


7 


yes 


29(7) 


19.8 




K+K+ 


35 


7 


yes 


16(8) 


5.0 




PP 


26 


6 


yes 


20(6) 


14.6 


p-Pb 


7T + 7T + 


21 


5 


yes 


29(9) 


24.8 




K+K+ 


26 


8 


yes 


29(9) 


23.1 




PP 


26 


8 


no 


20(8) 


7.6 



TABLE II. Parameters used in the reconstruction of the NA44 sources. The numbers in parentheses in the number of data 
points column is our estimate of the number of points which contain usable information. 



Ref. A R G [fm] X 2 /NDF 



K+K+ (S-Pb) 


7 


1 


0.92L0.08 


3.22L0.20 


53/31 


7T+7T+ (S-Pb) 


8 




0.46L0.04 


4.50L0.31 


18.1/16 


K+K+ (p-Pb) 


7 


1 


0.68L0.06 


1.71L0.17 


65/54 


TT+TV+ (p-Pb) 


8 


r 


0.38L0.03 


2.89L0.30 


16/25 



TABLE III. Gaussian fit parameters as obtained by NA44. 



Ref. A R D [fm] X 2 /NDF 



K+K+ (S-Pb) 


7 


1 


1.80L0.18 


2.64L0.22 


26/31 


7T+7T+ (S-Pb) 


8 




0.77L0.08 


3.54L0.33 


12.0/16 


K+K+ (p-Pb) 


7 


r 


1.10L0.10 


1.04L0.19 


55/54 



TABLE IV. Exponential fit parameters as obtained by NA44. 





{ PT } [MoV/c] 


A 


Rts [fm] 


Rto [fm] 


Rl [fm] 


7T+7T+ 


150 


0.56L0.02 


4.15L0.27 


4.02L0.14 


4.73±0.26 


7T+7T+ 


450 


0.55L0.02 


2.95L0.16 


2.97L0.16 


3.09±0.19 


K+K+ 


240 


0.82L0.14 


2.55L0.20 


2.77L0.12 


3.02±0.20 



TABLE V. Three dimensional Gaussian fit parameters from the S-Pb reaction as obtained by NA44. 
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FIGURES 
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2 4 6 



r 

FIG. 1. Sample plots of N^ h degree b-splines. In all panels, the knots are marked by carets and the knots at r = are 
actually Nb + 1 regular knots piled together. 



1.25 




50 100 50 100 150 

q [MeV/c] q [MeV/c] 

FIG. 2. Model proton correlation corresponding to (a) Gaussian proton source function and to (b) Dipole proton source. 
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f — I 

m 




10 20 30 40 10 20 30 40 10 20 30 40 
r [fm] r [fm] r [fm] 

FIG. 3. Reconstructions of the Gaussian source with different numbers of coefficients in the source image. In both panels, 
the model sources are the solid curves and the reconstructed source is the curve with the error band. 
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(b) 11 coefficients 




mm 



(c) 14 coefficients 



10-' 

10 20 30 40 10 20 30 40 10 20 30 40 
r [fm] r [fm] r [fm] 

FIG. 4. Reconstructions of the dipole source with different numbers of coefficients in the source image. In both panels, the 
model sources are the solid curves and the reconstructed source is curve with the error band. 
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10 20 30 40 10 20 30 40 10 20 30 40 
r [fm] r [fm] r [fm] 

FIG. 6. In three panels, the true dipole source is the solid curve and the reconstructed source is given by the points with 
errors. The knots in both panels are represented by carets. 
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5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 45 
r [fm] r [fm] 



FIG. 7. In both panels, the true Gaussian source is the solid curve and the reconstructed sources is given by the points with 
errors. In panel (a), the knots are evenly spaced between the limits of the imaging region and in panel (b), the "optimal knots" 
are used. The knots in both panels are represented by carets. Note that the source from panel (c) in Fig. ^ is reproduced here 
as panel (a). 




5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 45 
r [fm] r [fm] 



FIG. 8. In both panels, the true dipole source is the solid curve and the reconstructed sources is given by the points with 
errors. In panel (a), the knots are evenly spaced between the limits of the imaging region and in panel (b), the "optimal knots" 
are used. The knots in both panels are represented by carets. Note that the source from panel (c) in Fig. ^ is reproduced here 
as panel (a). 
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m 10 




5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 45 
r [fm] r [fm] 

FIG. 9. In both panels, the true Gaussian source is the solid curve and the reconstructed sources are the points with errors. 
In both reconstructions, the source is constrained to have zero derivative at the origin. In panel (a), we use 7 coefficients and, 
in panel (b), 9 coefficients. The knots in both panels are represented by carets. 




^ 10- 5 =r 



10 15 20 25 30 35 40 
r [fm] 



10 15 20 25 30 35 40 45 
r [fm] 



FIG. 10. In both panels, the true dipole source is the solid curve and the reconstructed sources are the points with errors. 
In both reconstructions, the source is constrained to have zero derivative at the origin. In panel (a), we use 7 coefficients and, 
in panel (b), 9 coefficients. The knots in both panels are represented by carets. 
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Legend: 
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- (a) 7t+tt + , S-Pb 200 GeV/c - 
Coulomb corrected 


- (b) TT+TT+, p-Pb 450 GeV/c - 
Coulomb corrected 






(c) K+K+, S-Pb 200 GeV/c 
Coulomb corrected 


(d) K+K+, p-Pb 450 GeV/c 
Coulomb corrected 
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FIG. 11. Pion, kaon and proton correlations for the S-Pb and p-Pb reactions. The points and narrow error bars correspond 
to the experimentally measured correlations. The histogram and wide error bars correspond to the restored correlations from 
the imaging analysis. The pion, kaon and proton correlations are from Refs. |9j, Jn, and |29(], respectively. 
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(c) K + K + , S-Pb 200 GeV/c 




(d) K + K + , p-Pb 450 GeV/c E 
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FIG. 12. Sources imaged from the S-Pb and p-Pb reactions. Where applicable, we have also plotted the Gaussian and 
dipole-shaped sources corresponding to NA44's fits. 
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FIG. 14. String fragmentation in the Lund string model. The hatched areas indicate regions of non-vanishing color field 
as well as the space-time region swept out by the oscillating strings. Points 1 and 2 indicate the space-time point where the 
hadron of interest breaks off of the main string. 
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